###########################################################################
###########################################################################
############# SCRIPT TO PRODUCE FIGURES                     ###############
###########################################################################
###########################################################################
###########################################################################

# Authors: 
#   
#   Yu Zeng (Asst. Prof, School of Government, Peking University)
#   Junyan Jiang (Asst. Prof, Dept of Political Science, Columbia University)
#   Jie Li (PhD candidate, Department of East Asian Studies, University of Vienna)
#   Christian Goebel (University Professor of Modern China Studies, Department of East Asian Studies, University of Vienna)

# Version: 1.0
# Tested under:
## R version 3.6.1
## Environment: Windows 11 x64, Intel i9-9980HK CPU, 32GB RAM

rm(list=ls())

require(plyr)
require(reshape)
require(reshape2)
library(haven)
library(ggplot2)
library(igraph)
library(Cairo)
library(data.table)
library(lubridate)
library(foreign)
library(grid)
library(gridExtra)
library(scales)
library(dplyr)
library(maptools)
library(rgdal)
library(sp)
library(sf)
library(raster)
library(plyr)

theme_set(theme_bw() + theme(text=element_text(family="Times"),legend.position="bottom",legend.direction="horizontal",axis.title.x=element_text(vjust=-.5),axis.title.y=element_text(vjust=1),plot.title=element_text(face="bold",size=17,vjust=1.5),axis.text.x=element_text(size=12),strip.text=element_text(size=13)))

windowsFonts(Times=windowsFont("TT Times New Roman"))

setwd("C:/Users/Ricardo/OneDrive/Documents/Research/ERCResponsiveness/Workspace/submission/")
d <- read_dta("replication_data.dta")



### Figure 1: Policy Map ### 

setwd("C:/Users/Ricardo/OneDrive/Documents/Research/ERCResponsiveness/Workspace/submission/map")

chinamap <- readRDS("chinamap.rds")

policy <- subset(d, year==2017, select = c(cityid, pinyin, firstyearHOA, fPolicyHome)) # fpolicyHome: pro-HOA policy from 09-18
beijing <-list(110000,'beijing',2010,1)
tianjin <-list(120000,'tianjin',2008,1)
shanghai <-list(310000,'Shanghai',2008,1)
chongqing <-list(500000,'chongqing',2009,1)
policy <-rbind(policy[!duplicated(policy$cityid), ],beijing)
policy <-rbind(policy,shanghai)
policy <-rbind(policy,tianjin)
policy <-rbind(policy,chongqing)

policymap <- left_join(chinamap, policy, by=c("CITY_ID" = "cityid"))

theme_set(theme_grey() + theme(text=element_text(family="Times"),legend.position="bottom",legend.direction="horizontal"))
p1 <-ggplot() +
  theme(plot.title = element_text(hjust = 0.5))+
  geom_sf(data = fortify(policymap),color="black",aes(fill=as.factor(fPolicyHome))) + 
  coord_sf(crs = "+proj=aeqd +lat_0=37 +lon_0=104")+
  scale_fill_manual(name="",
                    values=c("#FFFFFF","#000001","#00000"),
                    breaks=c(0, 1, NA),
                    labels=c("No PMR", "PMR", "Not a city/No information"))


l9 <- shapefile("dashline.shp")
china_map=rgdal::readOGR("bou2_4p.shp")

x <- china_map@data # read admin info
china_map1 <- fortify(china_map) # convert map to data frame

xs <- data.frame(x,id=seq(0:924)-1) # total of 925 shapes (islands included)

china_map_data <- join(china_map1, xs, type = "full")

windowsFonts(Times=windowsFont("TT Times New Roman"))

p2<-ggplot()+
  geom_polygon(data=china_map_data,aes(x=long,y=lat,group=group),color="grey40",fill="white")+ 
  geom_line(data=l9,aes(x=long,y=lat,group=group),color="black",size=0.5)+ 
  coord_cartesian(xlim=c(105,125),ylim=c(3,25))+ 
  theme(
    aspect.ratio = 1.25, 
    axis.text = element_blank(),
    axis.ticks = element_blank(),
    axis.title = element_blank(),
    panel.grid = element_blank(),
    panel.background = element_rect(fill = '#e5e5e5'),
    panel.border = element_rect(fill=NA,color="grey20",linetype=1,size=0.8),
    plot.margin=unit(c(0,0,-0.5,-0.5),"cm"))

tiff("policymap.tiff", res = 300,  width = 1800, height = 1800)
vie <- viewport(width=0.22,height=0.17,x=0.2,y=0.3)
p1
print(p2,vp=vie)
dev.off()


############ Figure 2: Validation of Independent Variable ##############

setwd("C:/Users/Ricardo/OneDrive/Documents/Research/ERCResponsiveness/Workspace/submission/")
source("stat_binscatter.r")

d$logprotest.r<-log(d$realestateprotest+1)
d$HomeownersComplt <- d$totalComplt*d$homeownersCompltPct
d$logpetition.r<-log(d$HomeownersComplt+1) 

library(miceadds)
f<-lm.cluster(logprotest.r~logpetition.r+logpop+loggdp,cluster="cityid",data=d)
summary(f)
txt1<-expression(beta == 0.25~"***")

g1<-ggplot(d, aes(y=logprotest.r,x=logpetition.r)) +
  stat_binscatter(bins = 10,geom="pointrange")+
  annotate(geom="text",x=1.25,y=3.65,label=as.character(txt1), size=6, parse=T)+
  xlab("Log Homeowner Complaint (MBL)")+ylab("Log Homeowner Protest (Wickedonna)")

d$homeowner_mgt <- rowSums(d[c("baiduindex_homecommittee","baiduindex_homebylaw","baiduindex_hoacondition","baiduindex_hoarule","baiduindex_hoaguideline","baiduindex_homeconvenant",
                            "baiduindex_propertycompany","baiduindex_proeprtyfee","baiduindex_propertylaw","baiduindex_pmr","baiduindex_propertyservice","baiduindex_propertymgt")], na.rm = TRUE)

d$log_homeowner_mgt <- log(d$homeowner_mgt+1)
f<-lm.cluster(log_homeowner_mgt~logpetition.r+logpop+loggdp,cluster="cityid",data=d)
summary(f)
txt2<-expression(beta == 0.38~"***")

g2<-ggplot(d, aes(y=log_homeowner_mgt,x=logpetition.r)) +
  stat_binscatter(bins = 10,geom="pointrange")+
  annotate(geom="text",x=1,y=4.1,label=as.character(txt2), size=6, parse=T)+
  xlab("Log Homeowner Petition (MBL)")+ylab("Log Baidu Index for Property owner-related terms")

tiff("main_validation.tiff", res = 300, height=1250, width=3000)
grid.arrange(g1,g2,ncol=2)
dev.off()


######################### Appendix ###########################

## Figure A.1

policy1 <- subset(policy, fPolicyHome==1, select = c(cityid, firstyearHOA, fPolicyHome)) # select cities with policy only
policy1 <- dplyr::mutate(policy1, year=ifelse(is.na(firstyearHOA),2018,firstyearHOA)) # create Column year: for all the existing policies (from 2009 to 2018)
policy1 <- subset(policy1, year>=2009, select = c(cityid, fPolicyHome, year))

ggplot(data=policy1, aes(x=year)) + geom_bar() + 
  scale_x_continuous(breaks=seq(2008,2019,1)) +
  scale_y_continuous(breaks=seq(2,20,2)) +
  xlab("")+ylab("No. of Cities Adopting Reform")


## Figure A.2： Trend of homeowners' complaints

HoYComplaints <- aggregate(HomeownersComplt~year, data=d, sum, na.rm=TRUE)
ToYComplaints <- aggregate(totalComplt~year, data=d, sum, na.rm=TRUE)

HoYComplaints <- HoYComplaints[which(HoYComplaints$year<2019), ]
ToYComplaints <- ToYComplaints[which(ToYComplaints$year>2007), ]

HoYComplaints_pct <- HoYComplaints$HomeownersComplt/ToYComplaints$totalComplt
HoYComplaints_pct <- cbind(HoYComplaints_pct, HoYComplaints$year)

colnames(HoYComplaints_pct)[2] <- "year"
plotComplaints <- data.frame(HoYComplaints_pct)

ggplot(data=plotComplaints, aes(x=year, y=HoYComplaints_pct)) +
  geom_line()+
  geom_point()+
  scale_x_continuous(breaks = 2008:2017)+
  labs(x="Year", y = "% Homeowners' Complaints")



